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We consider two regimes where a trapped Bose gas behaves as a one-dimensional system. In 
the first one the Bose gas is microscopically described by 3D mean field theory, but the trap is so 
elongated that it behaves as a ID gas with respect to low frequency collective modes. In the second 
regime we assume that the ID gas is truly ID and that it is properly described by the Lieb-Liniger 
model. In both regimes we find the frequency of the lowest compressional mode by solving the 
hydrodynamic equations. This is done by making use of a method which allows to find analytical 
oo ■ or quasi-analytical solutions of these equations for a large class of models approaching very closely 

the actual equation of state of the Bose gas. We find an excellent agreement with the recent results 
of Menotti and Stringari obtained from a sum rule approach. 



<4-l 
O 

CO 



T3 

o 



- 1—1 

X 



PACS numbers : 03.75.Kk, 05.30.Jp , 47.35.-H, 67.40.Hf 



I. INTRODUCTION 



Among the many experimental tools used to study Bose condensates of trapped cold atoms [1] , the study of collective 
modes has played a quite important role. Indeed they are of high importance both experimentally and theoretically. 
On the experimental side they provide direct in situ informations on the system, which are free of the quantitative 
interpretation of expansion experiments (this is of particular importance for dense systems where the standard mean 
field approximation does not work). On the theoretical side the low energy collective modes are the elementary 
excitations and as such they play an essential role in the physical understanding of these systems. This is clear from 
j-^. . the fact that they reduce to phonons for homogeneous systems. 

The recent developments in this field have shown, among others, a growing interest in two different directions. 
, One of them is toward the study of systems with reduced dimensionality [2]. Indeed it is experimentally possible 
to produce anisotropic trapping potentials which are strong enough in some directions to freeze the corresponding 
degrees of freedom. In these cases the ultracold atoms are in the ground state corresponding to the motion in these 
directions, since they have not enough energy to reach the related excited states. This produces systems which are 
effectively one-dimensional (ID) or two-dimensional (2D). These systems are of high fundamental interest since their 
physics is expected to have quite specific properties, different from those of 3D systems, which are in particular quite 
relevant for condensed-matter physics. By comparison atomic gases offer particularly clean, controllable and simple 
I . situations, quite close to model systems. 

Another trend is toward the study of strongly interacting systems [2]. Indeed most of the experiments have been 
performed with gases where the interactions are weak enough to be properly described by mean field theory. We 
define strongly interacting systems as those for which this description is no longer valid. In 3D this corresponds to 
dense gases. For Bose gases the coupling to molecular states has been observed and studied very recently in these 
dense systems [3], and one can hope to obtain in this way molecular condensates. Another interest is to bridge the 
gap between dilute Bose condensates on one hand, and liquid superfluid 4 He on the other hand where the interactions 
are very strong. For Fermi gases it is also quite interesting to go to dense systems because the BCS-like transition, 
which is presently very actively looked for, is expected to have a much higher critical temperature in this regime. 

Very recently Menotti and Stringari [4] (MS) have dealt with the problem of the collective oscillations in a ID Bose 
gas at zero temperature, trapped in a very elongated harmonic potential well. They considered both the case of a 
weakly and of a strongly interacting Bose gas. The first range covers the high to intermediate density regime, where 
mean field theory is always valid at a microscopic level and the system is described by the Gross-Pitaevskii equation 
[5]. The high density limit corresponds to the case where the Thomas- Fermi approximation of this equation is valid, 
and the elongated gas has physically the shape of a '3D cigar'. The intermediate density regime corresponds to the 
case where, for the transverse directions, all the particles are in the gaussian wavefunction describing the ground 
state of the transverse harmonic potential and only longitudinal degrees of freedom are left. This specific situation 
is called the 'ID mean field' regime. For all this range the system behaves for the low frequency modes as a ID 
system, since the transverse collective degrees of freedom correspond to much higher frequencies. However for this 
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effective ID behaviour the system is no longer effectively mean field (except in the intermediate density regime), i.e. 
specifically the chemical potential is no longer linear in ID density m, although the system is still microscopically 
described by 3D mean field. This is because one has to average over the transverse directions, which MS did from the 
Gross-Pitaevskii equation to obtain the chemical potential /x(ni) as a function of n\. 

In the second range studied by Mcnotti and Stringari [4], the gas goes from intermediate density (where the 
interaction is still weak and the ID mean field regime applies) to low density where the interaction is effectively 
strong. For all this range the gas is physically in a ID situation, with ID particle density m. A famous example 
of a ID strongly interacting system is the Lieb-Liniger [6] model of hard core bosons, of mass m, interacting via a 
repulsive delta potential giS(z), which they solved exactly. The paradoxical situation of a gas strongly interacting 
in a dilute limit (the so-called Tonks-Girardeau limit [7]) is actually due to the behaviour of the kinetic energy. A 
characteristic kinetic energy is h 2 k 2 /m, where k is typically related to the interparticle distance d by k <~ \ jd = n\. 
This is small in the dilute regime compared to a typical interaction energy n\gi, the ratio between these two energies 
being essentially the Lieb-Liniger coupling constant 7 = mgi/h 2 ni. For a D dimensional space the ratio behaves as 

n D 2 ^ D , so it docs not depend on density for D = 2 and decrease with density for D = 3. Specifically MS studied 
the Lieb-Liniger model with /J-{ni) obtained numerically from the Lieb-Liniger solution. In both ranges Mcnotti and 
Stringari used a sum rule approach to obtain the lowest frequency mode as a function of frequency. For the three 
limiting cases of the 3D cigar, the ID mean field and the Tonks-Girardeau limit, they recovered the results which 
they obtained from the analytical solution of the hydrodynamic equations. 

In the present paper we will use everywhere the hydrodynamic equations as a starting point. This is a quite natural 
and general approach since hydrodynamics is expected to be generically valid to describe the dynamics of the system 
for low frequencies and long wavelengths. In particular this approach is quite natural for dilute Bose condensates 
since hydrodynamics appear as a consequence of the Gross-Pitaevskii equation. Actually this is quite clear in 3D, 
but this link is more ambiguous in ID. Anyway we will assume that in ID hydrodynamics is a valid starting point to 
obtain the low frequency collective modes of the system. 

Quite recently two of us [8] (CL) considered how the fruitful hydrodynamic approach for collective oscillations in 
trapped Bose condensates, devcloppcd in the mean field regime [9], could be extended for strongly interacting systems 
where mean field is no longer valid. It was shown that, in this general case, the linear hydrodynamics could still 
be written in a quite convenient way. This made possible to find a number of specific functional dependence for 
/i(n) for which an analytical or quasi-analytical solution could be obtained. Conversely in the general case it appears 
possible to approximate [i(n) closely enough by some of these specific cases, considered as models, to obtain a very 
good approximation for the actual result, both for the frequency and the shape of the modes. This method can be 
applied to any mode. The flexibility of this modeling allows even in fact to invert experimental data covering a range 
of density to obtain the corresponding chemical potential n(n). Since this method is equivalent to solve in a simple 
and efficient way the hydrodynamic equations, it is the purpose of the present paper to apply it to the case of the ID 
Bose gas, in order to obtain the mode frequency from the solution of the hydrodynamic equations, rather than from 
the sum rule method that MS used. 

The paper is organized as follows. In the following section we will recall the CL method. Actually their original 
paper considered only explicitely the isotropic 3D case. So we will write the generalization to any dimension D, since 
we are interested in ID, and present the models we will use. Then in the next section we will consider the first case 
investigated by MS where the gas is microscopically 3D mean field, but behaves as a ID system for the low frequency 
modes. We will show that, even if mean field does not apply at a microscopic level, the same reduction to an effective 
ID problem can be obtained, by generalizing the procedure used by Stringari [10] in the case of mean field. We will 
obtain explicitely the ID effective chemical potential from the 3D n(n). Then, before turning to the case of the ID 
Bose gas, we will reinvestigate the case of the 3D mean field Fermi gas in an isotropic trap, which CL had already 
considered. This has been motivated by our preliminary results on the ID Bose gas which, although already quite 
reasonable, were not as good as expected. This led us to improve our method, essentially by making a first order 
correction to the result to take into account the small difference between the model ^i(n) and its actual value. In 
order to assess the results we have considered the 3D mean field Fermi gas, for which we had direct results from the 
numerical integration of the hydrodynamic equations. We have also checked the sum rule method in this case. The 
final result of all this is that we obtain the mode frequency with a relative precision of at least 10~ 3 and often much 
better. Naturally this precision, which was not looked for, is much better than necessary to compare with experiment, 
and it is likely that there will be some time before it proves useful to be so precise. On the other hand there is no 
reason to put aside this precision since we have it fairly easily. Moreover it gives us a very high confidence in our 
results, all the more because we can use different method which agree within this precision. Finally we will turn to 
the ID Bose gas, present our results and compare them with those of MS. The excellent agreement we find is a very 
good check for these two quite different procedures. 



2 



II. SOLVABLE MODELS FOR HYDRODYNAMICS 



The CL approach [8] can be extended quite generally to anisotropic traps. We will not dwell here on this general 
situation since we do not intend to make a specific use of it. We are only interested in very anisotropic traps. 
The intermediate situation of moderately anisotropic traps is indeed not very convenient theoretically to extract 
informations on the system. Experimentally it is also not so frequently used and one rather deals much more often 
with very elongated or very flat traps, which are quasi 1-dimcnsional or quasi 2-dimensional. Let us just mention 
that, in the a — p modeling (see below), we are restricted to a — 2 in this general situation. On the other hand it is 
possible to extend the use of quasipolynomial models (see below) to general anisotropic traps. We will rather show 
cxplicitely how the CL procedure can be used in any dimension D, having naturally in mind the interesting cases D=l 
and D=2. This will at the same time allow us to provide a short recall of this approach. 

As MS and CL we restrict ourselves to the case of the reactive hydrodynamics where dissipation is negligible and 
thermal effects can be omitted. This is valid for example at low enough temperature. In this situation we have just to 
write Euler equation mdv/dt = —Vjl(r,t) together with particle conservation dn/dt + V.(nv) = for density n(r,t). 
The global chemical potential jl(r,t) = jj,(n(r,t)) + V(r) has a contribution from the trapping potential V(r) and a 
contribution /i(n(r, t)) from the fluid itself, where /x(n) is the equilibrium dependence of the chemical potential on the 
density, as it results from thermodynamics. At equilibrium the particle density n (r) satisfies /i(no(r)) + V(r) = /i 
where fl is the constant value over the system of the global chemical potential. In particular jl is equal to the value 
of the trapping potential at the surface of the cloud (we take /i(n = 0) = 0) . For small fluctuations we introduce 
the departure of the chemical potential from its equilibrium value w(r, t) — jl(r,t) — fi = fi(n(r,t)) — /i(n (r)) = 
(dfi/dno)Sn(r 7 t) where Sn(r, t) = n(r,t) — n (r) is the density fluctuation. Assuming that the small fluctuations 
occur at frequency ui gives iivdn = V.(n v) = n V.v + v.Vn which, together with Euler equation iivmv = Vw, leads 
to: 

n V 2 w + Viio-Viu + muj 2 (dno/dfi) w = (1) 

Let us now restrict ourselves to an isotropic trap V(r) in D dimensions, where r is the distance from the origin. 
The equilibrium relation /z(no(r)) + V(r) = jl gives (d/j,/dno)(dno/dr) = —V'{r) where V'(r) is the derivative of V(r) 
with respect to r. Then for a mode with spherical symmetry w(r) = w(r) Eq.(l) becomes: 

2 

TTtU] T 

rw" + \D-l + rL'(r)]w' — L'(r)w = (2) 

V (r) 

where we have defined L(r) = ln(n (r)) with L'(r) = dL/dr. More generally a mode with an angular dependence 
w(r) = Yi m (9,ip)w(r) in 3D leads to: 

rw" + [D - 1 + rL'(r)]w' - + ^p-L'(r)]w = (3) 

r V'(r) 

where g(l) = 1(1 + 1). Similarly in 2D an angular dependence w(r) — exp(±ilip)w(r) gives the same equation with 
g(l) = I 2 . Finally in ID one gets again the same equation with g(l) = 0. In this last case it is better to rename the 
variable z since it is the abscissa and goes from —RtoR, where R is the cloud radius. It is then convenient to make 
explicit the dominant dependence at small r by setting w(r) = r l v(r) which leads to: 
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rv" + [2l + D-l + rL'(r)]v' - [y^ - l]L'{r)v = (4) 

This equation is valid in 3D and in 2D. It is also true in ID with I taking only two values : I = corresponding to 
the modes w(z) even with respect to z, or I — 1 corresponding to the modes odd with respect to z. 

We focus now on the most common case of the harmonic trap V(r) = ^mil 2 r 2 where Eq.(4) reduces to: 

rv" + [21 + D - 1 + rL'(r)]v' - (v 2 - l)L'(r)v = (5) 

with v 2 — uj 2 /il 2 . This equation has the quite convenient property to be scale invariant, if the same change of scale 
is naturally made for L(r). In particular it is unchanged under the replacement r — > r/R so we can take the cloud 
radius R as unity in the following. More generally the change of variable y = r a leads only to a modification of the 
constants in Eq.(5), provided again that the same change is made for L(r). One finds cxplicitely: 

d 2 v ,. dL.dv v 2 — ldL „ . . 

y d^ + iA + y ^d-y-—dy- V = ° (6) 
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where A = 1 + 2 ' + f ~ 2 ■ Finally we note that V does not change if the density n(r) is multiplied by a constant. So the 
absolute scale in density disappears and we have to deal only with the reduced density n(r) = n(r)/n(0). Similarly we 
introduce the normalized local chemical potential Jl(r) = /x(n(r))//x(n(0)) where /x(n(0)) is related to the gas radius 
R by /^(n(0)) = ^mVl 2 R 2 . This leads to Ji = 1 — r 2 , when r is expressed in units of R. 

When one takes the model dL/dy = — p/(l — y), Eq.(6) reduces to the hypergeometric differential equation: 

d 2 v dv v 2 — I 

The general solution of this equation, regular for y = 0, is the hypergeometric function F(a 1 b; A;y), with a and b 

defined bya + b = p + A — 1 and ab — ~P ]L —^- The boundary condition [8] at the surface of the cloud y = 1 give the 
further condition a = — n where n is a non negative integer. In this case the solution is just a polynomial. This leads 
for the normal mode frequencies to the explicit result: 

— =l + -n(n + p + A-l) (8) 

In particular for the lowest frequency mode which we will consider below, corresponding to I = and n = 1, the 
solution is merely v = (p + A)y — A, for frequency v 2 = a(l + A/p). The above model, which we call the a — p model 
in the following, corresponds explicitely to the equilibrium density n(r) = (1 — r Q ) p , which arises from the equation 
p, = 1 — (1 — n 1 /p) 2 / a linking chemical potential to density. In particular for a = 2, this is merely the power law 
jj = Actually, if we consider ID bosons, the 3D cigar, the ID mean field and the Tonks-Girardeau limit satisfy 

precisely this functional dependence [4] /2 = n 1 / p between chemical potential and density, with the respective values 
p = 2, p = 1 and p = 1/2. This leads to the following result for the even mode frequencies (/ = 0): 

^=2n+^(2n-l) (9) 



and: 



uj 2 n 



{)2 =(l + -)(2n + l) (10) 

for the odd mode frequencies (1 = 1). In this last case n = gives the 'dipole mode' (oscillation of the gas as a whole) 
with frequency u = CI. These results are in agreement with MS (their result gives the even or odd frequency modes, 
depending on the parity of their integer k). 

In addition to the simple a — p model considered above, CL found also for v(y) a wide class of quasi-polynomial 
solutions for more complex models. These quasi-polynomials are actually very rapidly converging series, which in 
pratice behave as polynomials (one can safely cut them off above some order) because the variable y is restricted by 
— V — 1- The corresponding models are dL/dy = — ^2^ =0 Pky k / (1 — y) arid can be considered as series expansion 
for (1 — y)dL/dy around y = with increased accuracy. They have in general K + 1 parameters in addition to a. In 
the following we will only use the K = 1 model (K = is the a — p model) which has three parameters a,po and p\ 
(we will call it the 3 parameters quasi-polynomial model). It is solution of: 



dv 

'^y2 ' ' i 1 " ' fiy 



with q 2 = -pi,qi = -(A +p ),9o = A,n = P\ ]L - zL and r = Po iL ^ i - The solution v = ^ n=0 a n y n is found by 
solving the recursion relation [8] between the a„'s resulting from Eq.(ll). Furthermore one requires a^v+i = which 
provides the equation allowing to find the mode frequency v 2 . The cut-off order N is chosen large enough to insure 
perfect convergence. In practice we have taken in our calculation N = 9 or 10, which is large enough for excellent 
convergence and small enough for very easy numerical calculations. 



III. REDUCTION FROM A 3D TO AN EFFECTIVE ID PROBLEM 



As we already discussed in the introduction there is one clear limit where a gas cloud can be considered as a one- 
dimensional system. This is the case where the transverse trapping potential is so strong, compared to temperature 
or interaction, that only degrees of freedom corresponding to motion along the weak longitudinal trapping direction 
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are left. The transverse degrees of freedom are completely frozen, the atoms being all in the ground state with respect 
to transverse motion. However there are weaker conditions under which, with respect to the modes, one has still an 
effective one-dimensional physics. Indeed if we deal with low enough frequencies, the transverse degrees of freedom 
will not be excited and we still have to deal only with the longitudinal degrees of freedom, which corresponds to a 
one-dimensional physical situation. This occurs for very elongated cigar-shaped traps where the transverse trapping 
is not too strong. This corresponds to the 3D cigar regime studied by Menotti and Stringari [4]. In their study MS 
used mean field theory to obtain the effective relation between chemical potential /j, and one-dimensional density m 
to be used in their ID treatment. Here we want to show how this result can be extended beyond mean field, and how 
in the framework of our approach we can derive the onc-dimensional equation for the modes in the 3D cigar regime. 
We basically follow the procedure of Stringari [10] which amounts to integrate Eq.(l) over the transverse variables. 
We consider an harmonic cigar shaped anisotropic potential, very elongated along the z axis: 

V{r^z)= 1 -m{ulrl+ulz 2 ) (12) 

with u> z <C u>±- The gas has a transverse Thomas- Fermi radius R(z) which depends on the location on the z axis. Since 
on the border of the cloud the local chemical potential /z(no(r)) is zero, the equilibrium relation between potential 
and chemical potential gives: 

l -m(ulR^z) + ^z 2 ) = ±mu 2 Z 2 (13) 

where Z is the maximal extension of the cloud along the z axis. We integrate Eq.(l) over the transverse position rj_, 
at fixed z, inside the disk rj_ < R(z). The two first terms of Eq.(l) are V(n .Vw) = V^(n .V^w) + V z (n .V z w). 
When the first term is integrated, from the divergence theorem it gives 27rno.Vj_w to be evaluated for r± = R{z). 
This is zero since the density no is zero at the cloud surface. 

Now comes the fact that the transverse degrees of freedom are not excited. In this case we have locally equilibrium 
in the transverse direction. This implies that the chemical potential fluctuation w(r) docs not depend on rj_ and 
depends only on z, so we have only to deal with w(z). Hence Vj_u> = (at lowest order in (uj z /uj±) 2 ). Then in 
order to integrate V z (n .V z w) over the transverse variable, it is more convenient to use, instead of r±, the local 
equilibrium chemical potential /z(no(r)) as a variable. From the equilibrium condition /x(no(r)) + V(r) = jj, we have 
dfi = —muj 2 _r±dr± and Eq.(l) leads us to: 

[ / d[i n a ]V 2 z w + V z [ / d[i n ].V z w + mLU 2 w / dfx— — =0 (14) 
Jo Jo Jo ufj, 

Here the integration over \x goes from zero (corresponding to the cloud border) to u(0, 0, z) = (j,(z) = \muj 2 (Z 2 — z 2 ) 
which is the local chemical potential on the z axis. Since w depends only on z we have taken it out of the integral, 
together with its derivatives. We are left with J dfi dno/dfi = no(/z(z)) where no(n(z)) is the equilibrium density on 
the z axis. Let us introduce: 

L{z) =ln[ ^ VnoW] (15) 
Jo 

The argument of the logarithm is just the transverse average / <iuno(/z) ~ fdrj_no(r) of the three-dimensional density. 
We have V z expL(2;) = no((J,(z))dfj,(z) / dz with d(j,(z)/dz = —muj 2 z z = —dV(z)/dz where V(z) — \muj 2 z 2 = V(0, z) is 
just the trapping potential on the z axis. Actually for our derivation we do not need an harmonic dependence on z 
for the trapping potential. On the other hand it is necessary to have an harmonic dependence on rj_, for our change 
of variables. Taking all this into account Eq.(14) becomes: 

2 
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w " + L'(z)w'--^L'(z)w = (16) 

This is exactly the equation we have already obtained for ID situations. In these cases we would have defined 
L(z) = lnno(z). Wc see that we have just to replace the one-dimensional equilibrium density no(z) by the transverse 
average n\{z) = Jdr± n (r) = 2tt Jd/ino(/j)/mw} of the three-dimensional density. This sounds a physically very 
reasonable result. In the mean field case we have fi = gno which leads to / d/zno(u) = p?{z)j2g and Eq.(16) becomes 
cxplicitely for a harmonic potential: 

j(Z 2 - z 2 )w" - zw' + v 2 w = (17) 
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with v = u 2 ju 2 z . This is as expected in agreement with Stringari [10]. The above mean field relation n\[z) — 
7r /j, 2 (z) j 'rrvJ\_g with g = Airh 2 a/m is also in agreement with MS result /i = 2hu>j_(ani) 1 / 2 for the 3D cigar regime. 
This makes clear that our above treatment works only for this 3D cigar regime, since the transition to the ID regime 
requires to include quantum effects, which are beyond the Thomas-Fermi approximation we had to use, consistently 
with hydrodynamics. Note finally that we could obviously apply the same treatment to a 2D pancake geometry. 



IV. FIRST ORDER CORRECTION 



As we have already explained in the introduction, we have improved the modeling used by CL and presented above 
by including a first order correction to the straight model evaluation of the mode frequency. This correction will be 
used in the two following sections, but we present it independently in this section for clarity. So let us now consider 
the possibility of correcting to first order the small difference between our model evaluation of the mode frequency, 
obtained for the model L (r) (this can be the a— p model or the 3 parameters quasi- polynomial model), and the actual 
result corresponding to the true L(r). This can be done conveniently for example by converting Eq.(6) into a second 
order differential equation which is formally identical to a Schrodingcr equation (that is without first derivative) . This 
is obtained by the change of function ip(y) = v(y)(y A n (y)) 1 ^ 2 which leads to: 



dy 



2 +V(y)i; = (18) 



with the effective potential: 



Eq.(19) is a Schrodinger equation corresponding to zero energy (and h 2 /2m = 1). If we have a change 5L = L — L , 
this gives a corresponding change for V(y) and a variation of the energy which can be calculated by conventional first 
order perturbation theory. In order to keep the energy equal to zero, we have also to give a compensating variation 
5v 2 for the frequency, which is just the correction we are looking for. Explicitcly this gives: 

The 8L" term can be integrated by parts (the integrated term is zero for the range of parameters we are interested 
in). When one uses L' = n /n and writes tp(y) in terms of v(y), one gets finally: 

r)i; 2 f 1 f 1 i/ 2 - 1 

— / dy(-n' Q )y A - 1 v 2 = / dy 8L> n ^v\- -v - yv'] (21) 

a Jo Jo a 



V. THE 3D MEAN FIELD FERMI GAS 



For the reasons presented in the introduction we study here, on the 3D mean field Fermi gas, how to improve the 
results obtained [8] by CL, concentrating specifically on the lowest frequency compressional mode. The first point 
to consider is how to best approach the actual p,(n) by a model p. mo d(n). This has to be done by minimizing some 
estimator / of the difference between our model and the actual physical equation. Since we have been interested 
in precision evaluation, the choice of the estimator / of the difference between our model and the actual physical 
equation of state becomes a relevant one. We have considered mainly three estimators. The two first ones are the 
pretty obvious: 

I 2 = / dn[p,(n) - fl mo d(n)} 2 and I 2 = / dp[n(p,) - n mod (fi)} 2 (22) 

Jo Jo 

and there is no reason why one of them should be systematically better than the other one. Note that one way to 
reconcile these two estimators would be to introduce a third one where, instead of dealing with the 'horizontal' or 
'vertical' distance between the two curves, one would consider the distance between the two curves perpendicular 
to them. However, since this is somewhat more complicated and does not bring a decisive advantage, we have not 
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implemented this solution. On the other hand all these estimators do not introduce a specific weighting while one 
could wonder for example if the high density regions are not more important than the low density ones. One possible 
way to introduce such a weighting is to require that the result of the first order correction we perform is as small as 
possible. However this correction Eq.(21) contains the shape of the mode itself, which is naturally not known before 
we have calculated the estimator and found its minimum. A possible way out of this dilemna is to retain only the 
factor in the correction which does not contain the mode and take for example: 



Jo 



dy y^no \SL'\ 



(23) 



However this estimator proved in our tests to be no better and even somewhat worse than the simpler estimators 
Eq.(22). In the end this question of the estimator choice turned out to be a secondary one because, once the first order 
correction is made, either of the estimators Eq.(22) is good enough and gives excellent precision. In the same spirit, 
once an estimator is chosen the minimization process may be difficult when there are different minima in competition 
(a well known problem in spin glasses for example). However because the corrected results are anyway very good, 
there is no need to go to the precise minimum, and an approximate minimization is enough. 

For the 3D mean field Fermi gas, with reduced units, the chemical potential is given in terms of the density by: 



3n 2 / 3 - 2Xn 
3-2A 



(24) 



where A is the coupling constant which goes from A = for the free Fermi gas to A = 1 when the attractive interaction 
is strong enough to produce a collapse. It is given by A = 2fcp|a|/7T in terms of the negative scattering length a and 
the Fermi wavevector kp. Naturally the interest of this study is that we have the results of the direct numerical 
integration of the hydrodynamic differential equation, which gives us a benchmark for our modeling. Even this direct 
integration has a limited precision, which in our case is not extremely high because there is no reason to work for 
very high precision. So CL used the classical Runge-Kutta routine with standard precision. The resulting numerical 
noise can be estimated directly from inspection of their results for v 2 as a function of A and its first order difference. 
We have an absolute precision of 10 -3 for these results on v 2 . 
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FIG. 1. Reduced lowest compressional mode frequency v 2 = (w/fi) 2 for a Fermi gas within the mean field approximation 
as a function of the coupling constant A. Full line: exact result from the numerical integration of the hydrodynamic equation. 
Long- dashed line: sum rule result. Short-dashed line: zeroth order a — p model. Dotted line: zeroth order 3 parameters 
quasi-polynomial model. 

The results of our calculations for v 2 as a function of A are summarized in Fig. 1 . It is clear that all the approximate 
methods give quite reasonable agreement with direct integration. So we have plotted in Fig. 2, with a much magnified 
scale, the difference between our various approximate calculations and the direct integration. As found by CL the 
zeroth order a — p model gives already quite satisfactory results. Indeed the maximum deviation from the direct 
integration is found around A w 0.9. For A = 0.9 it gives v 2 = 2.59 compared to 2.505 from direct integration. It is 
easy to see why this region for A is more difficult for a. — p modeling. This can be understood from the behaviour of 
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the equilibrium density no(r) near the origin. Since p, = 1 — r 2 we have in general no(r) ss no(0) + 0(r 2 ). On the 
other hand at the collapse A = 1 we have dfi/dno = and near the origin no(r) ~ no(0) + 0(r). Near the collapse 
no(r) is, so to speak, switching from one behaviour to another. So it is clear that its analytical behaviour will be 
more complex and its modeling will be accordingly more difficult. 



0.04 




-0.01 




FIG. 2. Difference Av 2 of the reduced mode frequency between various model results and direct integration of the hydrody- 
namic equation. Open circles: zeroth order a — p model. Open up-triangles: sum rule method. Open squares: 3 parameters 
quasi-polynomial model with first estimator in Eq.(22). Open down-triangles: 3 parameters quasi-polynomial model with sec- 
ond estimator in Eq.(22). Filled circles: corrected a — p model. Filled squares: corrected 3 parameters quasi-polynomial model 
with first estimator in Eq.(22). Filled down-triangles: corrected 3 parameters quasi-polynomial model with second estimator 
in Eq.(22). 

Making use of the 3 parameters quasi-polynomial model makes a very important improvement as can seen in Fig. 
1 and 2. Indeed for A ^ 0.4 the results are within the noise of the direct integration, that is within 10 -3 from the 
exact result. Then for higher A it starts to deviate with a maximum deviation of 0.03 for A = 0.9, then it gets back 
essentially to the exact value at the collapse. 

For comparison we have also plotted on Fig. 1 and 2 the result from the sum rule method [9] , which is basically a 
variational method. We have adapted it in the following way to our present problem. We start with the expression 
given by MS for the monopole frequency [4] 2uj^ 2 = — d ln(r 2 ) /<if2 2 (here we do not use the cloud radius as unit 
length). This relation is more conveniently expressed in terms of the oscillator length l a = \Jfijm£l as: 



J M 



d ln(r 2 ) 
dlo/lo 



(25) 



where, as above, v\ t — lo 2 m /SI 2 and the derivative is taken at constant number of particles N. Now, we express (r 2 ) as 
a function of l a , a and the coupling constant A defined above. In order to calculate (r 2 ) = J d 3 rr 2 no(r), we make in 
this integral the change of variables t = kpjk* with 3-7T 2 rio(r) = k 3 F and k* = %/2\a\. Using the equilibrium equation 
(24), we get r 2 = l\k* 2 P(X,t), with P(X,t) = A 2 - t 2 - |(A 3 - t 3 ). After a straightforward calculation, we get: 



S) = -±-l?k**F 3 (\) 
N=(l) = ^-l e o k* 6 F (X) 



(26) 
(27) 



where F q (X) 



J? dtt 3 p('?+ 1 )/ 2 dP/dt. From these two equations, we have: 

dln(r 2 ) = {F! 2 /F 2 )dX + 10^ 
d In N = (Fq/Fq) dX + 6^ 



(28) 
(29) 
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where F' — dF q /d\. Since we take the derivative at constant N in Eq.(25), we can eliminate dX from Eq.(28,29) and 
using Eq.(25), this yields: 



Vm 5-3F Fi/F'F 2 [M) 

The calculation of v 2 M (i.e the monopole frequency) therefore simply amounts to calculate four integrals. As it can be 
seen on Fig. 1 and 2, the sum rule method does remarkably well for A & 0.4 where it gives results within our direct 
integration noise. Then for larger A it starts to deviate and gets worse and worse when one goes to the collapse, 
where the deviation reaches 0.18. This is clearly linked [11] to the fact that the density has a linear dependence on r 
at the collapse, as we have seen, while the sum rule method is best suited when there is a quadratic dependence at 
the center. In this respect one has to keep in mind that, even near the collapse for the trapped gas cloud as a whole, 
most of the gas has a density far from the one corresponding to the collapse of an homogeneous gas, because of the 
inhomogcncity caused by the trapping potential. 

When we apply first order perturbative correction, as we have described above, we obtain a quite remarkable 
improvement. This is seen in Fig. 2 (the results have not been plotted in Fig. 1 since they would not be distinguishable 
from the direct integration ). Indeed the results we get from quasipolynomial modeling agree over the whole range 
of values for A with those of direct integration, within the numerical noise of 10~ 3 . This is quite gratifying, but 
perhaps not completely unexpected, taking into account that we have three adjustable parameters, plus a correction. 
More surprising is that the corrected a — p model gives also results which are almost within numerical noise, though 
the quasipolynomial results are slightly better. Since without correction the quasipolynomial model is clearly much 
better, as it can be seen in Fig. 1, this excellent final agreement of the a — p model is perhaps partially coincidental, 
although we find the same feature in the next section for the ID Bose gas. 

It is clear that the excellent agreement within 10 -3 of the results from the quasipolynomial and the a — p model 
gives us very much confidence in the validity of the result itself within this precision. We could ignore the direct 
integration and only infer the final result from what we obtain from the various models. This is naturally the point 
of view that we will adopt in the next section when we will study ID bosons. 




FIG. 3. Lines of equal values in the a — p plane for the first of the estimators Eq.(22). Starting from the minimum which 
is at 3.1CT 5 , the black ellipse corresponds essentially to 3.1(T 5 ,the dashed lines to 4.1CT 5 , 5.1CT 5 , 6.KT 5 , 8.1CT 5 and 9.1CT 5 
respectively for I 2 . The band of open squares corresponds to the region where the reduced frequency calculated by the corrected 
a — p model is within ±1CP 3 from the direct integration. This band is also shown in Fig. 4. 



In this respect, as we mentionned already, it is interesting to note that the first order perturbative correction takes 
essentially care of the difference which appears when we use different modeling and/or use different estimators and/or 
use different minima of these estimators. In particular, as we have already mentionned, if one requires very precise 
results, it is not always easy to locate the absolute minimum of an estimator. The use of the first order perturbative 
correction solves this problem since, as we have seen, the dispersion of the corrected results is very small compared 
to the one of the uncorrected ones. So we are allowed a little imprecision in the choice of the parameters of our 
modeling, since the first order perturbative correction will compensate for the resulting error in the mode frequency. 
It is naturally interesting to study this point in more details. We have done it for a worst case situation A = 0.9, a 
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value near a maximal dispersion of the results as we have seen. We have only considered the a — p model in order 
to have only a two-dimensional parameter space. We have covered the range 1.15 < a < 1.4 , 1.65 < p < 1.9. The 
lines of equal values for the first of the estimators Eq.(22) are drawn in Fig. 3. Since we are near the minimum these 
lines arc elliptical. One sees that the minimum is quite shallow for the direction in the a — p plane corresponding to 
the large axis of these ellipses. We have also shown in this figure the region where the corrected a — p model gives a 
result within ±10 -3 from the direct integration. It is quite satisfactory that this region goes very near the minimum 
of our estimator (the difference is basically within the noise) and is essentially oriented along the large axes of the 
ellipses (naturally we have generically a line in the a — p plane corresponding to any given value of v 2 ). Again for 
the corrected a — p model, we have plotted in Fig. 4 the lines of equal values for v 2 being respectively at ±10 -3 , 
±2.10 -3 , ±5.10 -3 and ±7.10 -3 from the direct integration result v 2 = 2.505. We see that one can take values for 
(a,p) which correspond to an estimator three times larger than the minimum and still obtain a result for v 2 which is 
at worst only within ±5.10 from the exact result. The corresponding value of the estimator shows that the average 
deviation between the model and the exact /z(n) is roughly 1%. This shows that, even with a modeling which is not 
very precise, we can obtain a very good result for v 2 . 




FIG. 4. Lines of equal values in the a — p plane for v 2 being respectively at ±10 -3 (dotted line), ±2.1CF 3 (short dashed- line) , 
±5.10 _3 (long dashed-line) and ±7.10 -3 (full thick line) from the direct integration result v 2 — 2.505. The full thin line 
corresponds to values in the a — p plane where the first order correction to the a — p model is actually zero. 

When one sees the range of values for the parameters giving essentially the correct mode frequency, one may wonder 
if there is not another criterium than the above estimators to find a priori the best set of parameters. An interesting 
possibility in this direction is to say that the best parameters should not require any first order correction at all. So 
one should use sets of parameters giving zero correction. These sets correspond to a line in the a — p plane which is 
shown in Fig. 4. We have calculated the corresponding mode frequencies. We have found that, although the results 
are fairly independent of the parameter set, they fall systematically slightly above, at 2.514 , compared to the correct 
result 2.505 . Actually either of the estimators Eq.(22) turns out to be the best criterium to find the correct result 
(the difference between these two being within the numerical noise). This is seen in Fig. 2 where both estimators 
have been used. 



VI. THE ID BOSE GAS 

We apply now the previously developed method to study the hydrodynamic modes of the trapped ID Bose gas with 
repulsive interactions. At zero temperature, the interactions arc characterized by a single parameter: the 3D s-wave 
scattering length a. The system we consider is the same as the one studied by Menotti and Stringari [4]. Namely, 
this is a Bose gas in a very anisotropic trap (u)± S> lo z ) for various one-dimensional configurations. The transverse 
oscillator length a± — ^/h/muj± is supposed to be much larger than the scattering length a± ^> a. As we mentionncd 
already in the introduction, MS identified three limiting regimes corresponding physically to decreasing ID densities 
ni(0) at the center of the trap. At high density (m(0)a » 1), the system is in the radial Thomas-Fermi regime (3D 
cigar). In this case the chemical potential is related [4] to the ID density by \x = 2fiLj^{ani) 1 / 2 , as indicated already 
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in section III. In the second regime of intermediate densities, it is convenient to introduce the effective ID scattering 
length a\ related to the Lieb-Linigcr coupling constant g\ by g\ — 2h 2 jma\. In the limit a±_ 3> a, Olshanii [13] 
showed that a\ — Oj_/a [14]. In this intermediate density regime defined by m(0)ai > 1 > m(0)a, the gas is a ID 
quasi-condensate (ID mean-field regime) and the chemical potential is given by /i = 2huj±ani (with our conventions \x 
is zero for zero density). Finally at low density 1 » m(0)ai, the gas consists of impenetrable bosons (Tonks-Girardeau 
regime) where fi = -K 2 fi 2 n\/2m. 

Our needed input, in order to obtain the mode frequencies, is the equation of state /u(m) for the homogeneous gas 
at zero temperature. For the transition between the 3D cigar and the ID mean-field regimes (high density domain), 
it can be obtained from the numerical solution of the 3D Gross-Pitaevskii equation in a cylindrical geometry [4]. 
Actually, it appears that the numerical results [12] can be very well approximated analytically by: 

/u(m) = Twjjl{\/1 + 4ma - 1) (31) 

This formula gives the exact behavior in the two limiting cases (3D cigar and ID mean-field) defined above and 
interpolates very well in between (to better than 2%). In the following, we will use this approximate analytic 
expression which is easier to manipulate and avoids numerical problems encountered when the equation of state is 
only known for a fixed set of points. In the high density domain, we use the dimensionless variable n\(0)a to go from 
one regime to the other. 

The equation of state for the transition between the ID mean-field and Tonks-Girardeau regimes (low density 
domain) can be obtained from the Lieb-Liniger solution [6] . Lieb and Liniger gave a closed expression for the energy 
per particle as a function of the density in the form of an integral equation. It shows that the chemical potential is a 
universal function of nidi, which has to be evaluated numerically. Mcnotti and Stringari [4] calculated this equation 
of state and made their result available [12]. In the following, we use their data. In the low density domain, the 
dimensionless variable we use is m(0)ai. 

Using these equations of state, we compute the mode frequency using the following procedure: for each value of 
ni(0)a (or m(O)ai) we fit the equation of state (for m varying between and ni(0)) with an analytic model (either 
a — p or 3 parameters quasi-polynomial) using one of the estimator of equation (22); the zero order mode frequency 
is then obtained by inserting the value of the best set of parameters in the formula giving the mode frequency for the 
analytic model (Eq.(8) for the a — p model); we then compute the first order perturbation correction to the mode 
frequency. 

We first discuss the different models used to compute the (squared) reduced mode frequency v 2 = ui 2 /ui 2 . Fig. 5 
shows v 2 as a function of ni(0)a for a system where a±/a = 100. Four curves are plotted. They correspond to the 
a — p model, the quasi-polynomial model, and to the same two models corrected to first order of perturbation theory. 
Actually, the 3 parameters quasi-polynomial model is used with only two parameters by fixing a = 2. We checked 
that, in the case of the ID Bose gas and unlike that of the 3D Fermi gas, it does not make a significant difference to 
let the three parameters free or to set a = 2, whereas it is much faster and easier to work with only two parameters. 

V 2 '«A ' ' ' ' ' ' ' ' ' 1 
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2.5 

-6 -5-4-3-2-1 1 2 3 4 

log(n,(0)a) 

FIG. 5. Lowest compressional mode of the trapped ID Bose gas with a±/a — 100. The (squared) mode frequency v 2 — to 2 /lu 2 
is plotted as a function of log (m(0)a). The crosses (+) correspond to the a — p model; the stars (★) to the corrected a — p 
model; the circles (o) to the quasi-polynomial model; and the diamonds (o) to the corrected quasi-polynomial model. Actually 
at this scale all the symbols fall on top of each other, except for the crosses. 
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There are much more points plotted between the 3D cigar and the ID mean-field regimes (high density domain) 
than between the ID mean-field and the Tonks- Girardeau regime (low density domain). This is a mere consequence 
of the fact that we used an analytic expression for the equation of state in the first domain Eq.(31), whereas in 
the second one we used the numerical data of MS. Having a fixed set of points for the equation of state makes the 
numerical evaluation of integrals (needed to compute the correction 8v 2 , for example) more difficult because we are 
forced to use a primitive integration algorithm and, therefore, puts a limit on the precision and on the number of 
points that can be calculated safely. The fact that the junction between the two domains is not perfect (in the region 
of the ID mean-field regime) is a consequence of taking a finite value for a±/a instead of the limit a±/a —> oo. As 
noted by MS, the ID mean-field regime can only be identified provided aj_/a 3> 1. 

In order to compare our various calculations, we look now at our results with a much larger scale, which amounts 
to using a magnifying glass. We take as reference the results of the corrected quasi-polynomial model. As discussed 
in the preceding section we expect this reference to be the best of our results and to be very precise. The difference 
between this reference and the three first calculations (a—p, quasi-polynomial and corrected a — p models) are plotted 
in Fig. 6. Two points already made when discussing the 3D Fermi gas are worth stressing again in the case of the 
ID Bose gas. First, once corrected, the results of the a — p and the quasi-polynomial models agree remarkably well. 
Second, their agreement is at the absolute 10 -3 level for v 2 . 

The a — p model is exact in the three limiting regimes: a — 2 and p = 2 in the 3D cigar regime, which gives 
v 2 = 5/2; a — 2 and p = 1 in the ID mean-field regime, which gives v 2 = 3; and a = 2 and p = 1/2 in the 
Tonks-Girardeau regime, which gives v 2 = 4. In between these limits, it is apparent on Fig. 5 and 6 that the a — p 
model has some difficulties in predicting precisely the correct value of the mode frequency. This is what originally 
motivated the use of the 3 parameters quasi-polynomial model and the development of corrections using perturbation 
theory. This difficulty can actually be understood in much the same way as for the 3D Fermi gas. For example in the 

1 /2 

transition between the 3D cigar and the ID mean field, one finds that for large ni(0)a one has essentially p, ~ 

for most of the n\ range, while for small m this dependence turns into ~ n\. This switch of analytical behaviour 

is difficult to follow for the simple a — p model, hence the somewhat unsatisfactory result. 

For the a — p and quasi-polynomial models, using an analytic expression for the equation of state (high density 
domain) gives better results than using numerical data (low density regime), as can be seen on Fig. 6. This is linked 
to the difficulty, mentionned above, of using purely numerical data as an entry for However, the corrected 

models agree at the same level in both domains. 
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FIG. 6. Lowest compressional mode of the trapped ID Bose gas with a±/a — 100. The (squared) mode frequency v 2 — ui 2 /wl 
calculated using the corrected quasi-polynomial model is taken as a reference. The difference between this reference and three 
different calculations (a—p, corrected a— p and quasi-polynomial models) are plotted as a function of log (ni(0)o). The crosses 
(+) correspond to the a — p model; the stars (*) to the corrected a — p model; and the circles (o) to the quasi-polynomial 
model. 

We will now compare our results with those of Ref. [4], first discussing the high density domain. Menotti and 
Stringari calculated v 2 as a function of Naa±/a 2 , where a z — yjTi/mu) z is the longitudinal oscillator length and N the 
number of particles. In this paper, we calculate the same quantity as a function of the more convenient m(0)a, which 
is an increasing function of Naa^/a 2 . Finding the relation between these two numbers amounts to obtaining the 
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normalized density profile. This can be done analytically in the two limiting regimes [4] and numerically in between. 
This allows to plot the MS result for the mode frequency and our corrected quasi-polynomial result (which we expect 
to be the most precise of our results) on the same graph, see Fig. 7. It shows excellent agreement, at the 10~ 2 
level. Asking for agreement at a better level is not sensible here, as we used an analytic expression approximating the 
equation of state within a few percent and not the exact equation of state. This shows up in the fact that for high 
densities (ni(O)a > 1), the corrected quasi-polynomial seems to predict values of the mode frequency that are higher 
than the result of MS. This would be in contradiction with the fact that the sum rule approach they used is known 
to be an exact upper bound [9]. 



High density domain 
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FIG. 7. Lowest compressional mode of the trapped ID Bose gas in the high density domain. The (squared) mode frequency 
v 2 = Ld 2 /io 2 is plotted as a function of log(ni(0)a). The diamonds (o) correspond to the corrected quasi-polynomial model. 
The full line is the result of Menotti and Stringari. 

As a further check on our method, we also calculated the expected frequency for two experiments. For this purpose, 
the exact equation of state was used, in the form of the numerical solution of the 3D Gross-Pitaevskii equation [12]. 
For the experiment of Ref. [15], where m(0)o = 0.22, we hnd v 2 = 2.845 ± 0.003 and v 2 = 2.911 ± 0.002 for the 
experiment of Ref. [16], where ni(0)a = 0.10. This is in complete agreement with the values first obtained by MS. 



Low density domain 
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FIG. 8. Lowest compressional mode of the trapped ID Bose gas in the low density domain. The (squared) mode frequency 
v 2 = lo 2 /lj 2 is plotted as a function of log (m(0)ai). The diamonds (o) correspond to the corrected quasi-polynomial model. 
The full line is the result of Menotti and Stringari. 

We consider finally the low density domain. As in the high density domain, we first have to relate the dimensionlcss 
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variable we use ni(0)ai to Na\/a 2 z used in Ref. [4]. This is done analytically in the ID mean-field and in the Tonks- 
Girardeau regime [4] and numerically in between. The mode frequency obtained by the corrected quasi-polynomial 
model and that calculated by MS are plotted in Fig. 8. The agreement is again very good. 

The results of MS are in very good agreement with our result over the whole range of n\{Q)a. In the case of the ID 
Bose gas, the sum rule approach [9] seems to give not only an upper bound but to be quite near the exact result. As 
already discussed in the case of 3D Fermi gas, this is related to the fact that the density has a quadratic dependence 
on r near the center of the trap. In other words, the parameter a is always very close to 2. 



VII. CONCLUSION 



In this paper we have calculated, as a function of density at the center of the trap, the frequency of the lowest 
compressional mode of a ID trapped Bose gas, taking the reactive hydrodynamic equations as a starting point. We 
have considered two density regimes. In the first one the density decreases from high to intermediate, but the Bose gas 
is always described microscopically by 3D mean field theory. In the high density region the Bose gas has the shape of 
a 3D cigar, while in the intermediate density region the atoms are in the one-particle ground state with respect to the 
transverse motion and the gas behaves as a ID system. Nevertheless for all this regime the low frequency modes are 
described by ID hydrodynamic equations, but the effective equation of state Eq.(31) is no longer given by mean field. 
The second regime that we have considered goes from intermediate to low density and we have taken the Lieb-Linigcr 
model to describe it, which corresponds to a ID Bose gas going from a weakly to a strongly interacting situation. 

In order to solve the hydrodynamic equations we have made use of a very recent approach which allows to find 
analytical or quasi-analytical solutions of these equations for a large class of models. These model solutions allow to 
approximate very nearly the exact equation of state which is the only input of the hydrodynamic equations coming 
from the physical properties of the system. When in addition a first order correction has been made, we have been 
able to check that this method gives the correct mode frequency with at least a relative precision of 10~ 3 which is 
more than necessary for any practical purpose. On the other hand the simplest of this model gives a very easy and 
convenient analytical solution. Taken together the ensemble of these models allow to cover all the range from simple 
analytical solutions to very precise numerical solutions. We have compared our results to those obtained by Menotti 
and Stringari from a sum rule approach, and we have found an excellent agreement. 

The method used in this paper is quite powerful. It is not restricted to the lowest frequency mode and it can be 
used as well for any higher frequency mode. We have not presented the corresponding results in the present paper 
only to avoid to make it oversized, but this would have been quite easy. Another interest of our method is that it is 
not restricted to mean field and can be applied to any equation of state corresponding possibly to a very dense Bose 
gas. For example we could very well consider the situation where the gas is dense enough so that the ID intermediate 
density situation can no longer be described by mean field [17,18]. Moreover the convenience of analytical solutions 
makes it possible to invert the method [8] and to extract the effective equation of state of the gas from experimental 
data on the variation of the modes frequencies as a function of the density. 
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